

capture cd "~/Dropbox/Projects/Zimbabwe/Protected Villages and Political Outcomes/Replication"
use "data/Zimbabwe AB all coding.dta", clear

quietly {
	
	merge m:1 zimwardsid using "data/zim_geo_wards.dta"

	drop if _merge == 2
	drop _merge

	replace ward_centroid_x = x
	replace ward_centroid_y = y
}

****************************
***   Create variables   ***
****************************


quietly {

	egen round_group = group(round)
	*gen log_area = log(wardareasq)
	gen pv = 0
	replace pv = 1 if pv_area != ""
	replace pv = 1 if pv_border == 1
	encode pv_area, gen(pv_code)
	
	gen rural = strpos(local_auth, "RDC") > 0
	replace rural = . if local_auth == ""

	gen log_moz_dist = log(dist_to_moz)
	gen log_harare_dist = log(dist_to_harare)
	gen log_provcap_dist = log(dist_to_provcap)

	
	* Treatment variables for diff in disc
	gen dist_to_border = native_dis
	replace dist_to_border = -native_dis if native == 0
	gen treatment = (pv == 1 & dist_to_border > 0) if pv != . & dist_to_border != .
	gen treatment_int = treatment * dist_to_border
	gen native_int = native * dist_to_border
	gen pv_int = pv * dist_to_border
	label var treatment "Effect of PV"
	
	* For cohort effects
	gen post = 0 if year_birth != . & year_birth >= 1945
	replace post = 1 if year_birth >= 1980 & year_birth != .

	* Create indices
	alpha good_perf_pres good_perf_mps, std gen(gov_performance)
	alpha sat_manag_econ sat_manag_poverty sat_creating_jobs sat_prices_down sat_narrowing_inc_gap sat_ensuring_eat, std gen(sat_econ)
	alpha sat_reducing_crime sat_impr_health_ss sat_address_educ_needs sat_prov_water_and_sanit_ss sat_fight_corrupt sat_maintain_roads sat_prov_electricity, std gen(sat_issues)
	alpha reject_*, std gen(reject_auth)
	alpha close_inc_party trust_inc trust_pres, std gen(trust_ruling)
	alpha trust_police trust_army, std gen(trust_security)
	alpha trust_courts trust_ec, std gen(trust_inst)
	alpha att_rally-work_for_cand registered voted, std gen(formal_politics)
	alpha radio-social_media, std gen(consume_politics)
	alpha resp_suspicious resp_hostile nbr_suspicious nbr_fear, std gen(fear)
	alpha reps_approached enum_threatened enum_physical, std gen(enum_threat)
	alpha post_office school police_station health_clinic pavement, std gen(infrastructure)
	alpha electricity piped_water sewage, std gen(services)
	
	label var gov_performance "Gov performance"
	label var sat_econ "Satisfied (econ)"
	label var sat_issues "Satisfied (issues)"
	label var reject_auth "Reject auth. regimes"
	label var trust_ruling "Trust ZANU-PF"
	label var trust_security "Trust security"
	label var trust_inst "Trust institutions"
	label var formal_politics "Formal politics"
	label var consume_politics "Consume politics"
	label var fear "Resp/neighbor fear"
	label var enum_threat "Enum. threatened"
	label var infrastructure "Infrastructure in EA"
	label var services "Services in EA"

}

*****************************************************************************************

*** MAIN TEXT

*****************************************************************************************

* Define outcomes

global govperform gov_performance sat_econ sat_issues
global economic no_job infrastructure services
global dem_support reject_auth parties_needed two_term democracy
global zanu trust_ruling trust_security trust_inst
global polactivity formal_politics att_demo disc_pol consume_politics
global community assoc_comm att_comm_meet raise_issue should_agreecommunity
global state_capacity fear enum_threat police_vehicles
global dep_outlook close_main_opp_party trust_opp dem_doesntmatter


************************
***   Main results   ***
************************

set scheme plotplainblind

foreach macro in govperform economic dem_support zanu {

	capture mkdir "plots/`macro'"
	
	foreach var of varlist $`macro' {
	
	est clear
	loc j=0
	
	disp "`var'"
	
		forvalues i = 5(0.5)9.5 {
		
		disp "band = `i'"
		quietly {

			eststo: areg `var' c.treatment native dist_to_border pv ///
			treatment_int native_int pv_int i.round_group year_birth ///
			if dist_to_moz < 150 & abs(dist_to_border) < `i', ///
			absorb(province) cluster(local_auth)
			
			}
			
			est sto m`j'
			loc `++j'
			
		}

		coefplot (m0 m1 m2 m3 m4 m5 m7 m6 m7 m8 m9), pstyle(gs15) mcol(black) ///
		keep(treatment) asequation ///
		yline(0, lcolor(black) lwidth(0.25)) levels(90 95) msize(2) symbol(circle) ///
		ciopts(lwidth(0.8 0.5)) yscale(range(-5, 5)) ///
		ylabel(-5(1)5,labcolor(black) labsize(large)) ///
		xlabel(,labcolor(black) labsize(large)) ///
		ytitle("Effect size", axis(1) color(black) size(large)) ///
		xtitle("Bandwidth (km)",  color(black) size(large)) legend(off) ///
		coeflabels(, wrap(15)) ///
		swapnames eqrename (m0 = "5" m1 = "5.5" m2 = "6" m3 = "6.5" m4 = "7"  ///
		m5 = "7.5" m6 = "8" m7 = "8.5" m8 = "9" m9 = "9.5") vertical

		quietly graph display Graph, ysize(1) xsize(1)
		graph export "plots/`macro'/`var'.pdf", replace

	}	
	
}


************************************
***   Alternative explanations   ***
************************************

set scheme plotplainblind

foreach macro in polactivity community state_capacity dep_outlook {

	capture mkdir "plots/`macro'"
	
	foreach var of varlist $`macro' {
	
	est clear
	loc j=0
	
	disp "`var'"
	
		forvalues i = 5(0.5)9.5 {
		
		disp "band = `i'"
		quietly {

			eststo: areg `var' c.treatment native dist_to_border pv ///
			treatment_int native_int pv_int i.round_group year_birth ///
			if dist_to_moz < 150 & abs(dist_to_border) < `i', ///
			absorb(province) cluster(local_auth)
			
			}
			
			est sto m`j'
			loc `++j'
			
		}

		coefplot (m0 m1 m2 m3 m4 m5 m7 m6 m7 m8 m9), pstyle(gs15) mcol(black) ///
		keep(treatment) asequation ///
		yline(0, lcolor(black) lwidth(0.25)) levels(90 95) msize(2) symbol(circle) ///
		ciopts(lwidth(0.8 0.5)) yscale(range(-5, 5)) ///
		ylabel(-5(1)5,labcolor(black) labsize(large)) ///
		xlabel(,labcolor(black) labsize(large)) ///
		ytitle("Effect size", axis(1) color(black) size(large)) ///
		xtitle("Bandwidth (km)",  color(black) size(large)) legend(off) ///
		coeflabels(, wrap(15)) ///
		swapnames eqrename (m0 = "5" m1 = "5.5" m2 = "6" m3 = "6.5" m4 = "7" m5 = "7.5" ///
		m6 = "8" m7 = "8.5" m8 = "9" m9 = "9.5") vertical

		quietly graph display Graph, ysize(1) xsize(1)
		graph export "plots/`macro'/`var'.pdf", replace

	}	
	
}



*****************************************************************************************

*** APPENDIX

*****************************************************************************************

******************************
***   Summary statistics   ***
******************************

capture mkdir "tables"

est clear


global govperform gov_performance sat_econ sat_issues
global economic no_job infrastructure services
global dem_support reject_auth parties_needed two_term democracy
global zanu trust_ruling trust_security trust_inst
global polactivity formal_politics att_demo disc_pol consume_politics
global community assoc_comm att_comm_meet
global state_capacity fear enum_threat police_vehicles
global dep_outlook close_main_opp_party trust_opp dem_doesntmatter


preserve
keep if dist_to_moz < 150

estpost summ native pv $zanu $govperform $dem_support $polactivity $community $economic $state_capacity $outlook

esttab using "tables/sumstats.tex", replace cells("count mean(fmt(3)) sd(fmt(3)) min max") label nomtitle noobs nonumber
restore


***********************
***   Full tables   ***
***********************

capture erase "tables/tables.tex"


foreach macro in zanu govperform dem_support polactivity community economic  state_capacity dep_outlook {
	
	est clear
		
	quietly: foreach var of varlist $`macro' {
	
		disp "`var'"
		
		* MSE
		
		rdbwselect `var' dist_to_border if dist_to_moz < 150 & pv != ., ///
			bwselect(msetwo)
		local bandleft=e(h_msetwo_l)
		local bandright=e(h_msetwo_r)
		
		
		eststo: areg `var' treatment native dist_to_border pv treatment_int ///
			native_int pv_int i.round_group year_birth if dist_to_moz < 150 & ///
			dist_to_border <= `bandright' & dist_to_border >= -`bandleft', ///
			absorb(province) cluster(local_auth)
			
	

	}
	
	esttab ///
	using "tables/tables.tex", b(3) se(3) ar(3) ///
	stats(N, fmt(0 3 3 3 3 "")  ///
	labels(`"Observations"' `F-stat')) ///
	star(* 0.1 ** 0.05 *** 0.01) ///
	mtitles(`var') ///
	keep(treatment) title("`macro'" \label{tab:`macro'}) style(tex) booktabs label append

			
}



*******************************
***   Geography (balance)   ***
*******************************

set scheme plotplainblind

global land_type share_purchase share_ttl share_undesignated share_urban share_resettle share_crop
global geographic share_water2 share_lake share_forest share_natpark share_shrub share_grassland share_sp_vegetation share_wetland share_bare 
global distance log_moz_dist log_harare_dist log_provcap_dist

foreach macro in land_type geographic distance {

	capture mkdir "plots/`macro'"
	
	foreach var of varlist $`macro' {
	
	est clear
	loc j=0
	
	disp "`var'"
	
		forvalues i = 5(0.5)9.5 {
		
		disp "band = `i'"
		quietly {

			eststo: areg `var' c.treatment native dist_to_border pv ///
			treatment_int native_int pv_int i.round_group year_birth ///
			if dist_to_moz < 150 & abs(dist_to_border) < `i', ///
			absorb(province) cluster(local_auth)
			
			}
			
			est sto m`j'
			loc `++j'
			
			*tab native pv if dist_to_moz < 150 & abs(dist_to_border) < `i'
			*noisily lincom treatment
		}

		coefplot (m0 m1 m2 m3 m4 m5 m7 m6 m7 m8 m9), pstyle(gs15) mcol(black) ///
		keep(treatment) asequation ///
		yline(0, lcolor(black) lwidth(0.25)) levels(90 95) msize(2) symbol(circle) ///
		ciopts(lwidth(0.8 0.5)) yscale(range(-5, 5)) ///
		ylabel(-5(1)5,labcolor(black) labsize(large)) ///
		xlabel(,labcolor(black) labsize(large)) ///
		ytitle("Effect size", axis(1) color(black) size(large)) ///
		xtitle("Bandwidth (km)",  color(black) size(large)) legend(off) ///
		coeflabels(, wrap(15)) ///
		swapnames eqrename (m0 = "5" m1 = "5.5" m2 = "6" m3 = "6.5" m4 = "7" m5 = "7.5" ///
		m6 = "8" m7 = "8.5" m8 = "9" m9 = "9.5") vertical

		graph display Graph, ysize(1) xsize(1)
		graph export "plots/`macro'/`var'.pdf", replace

	}	
	
}


****************************************
***   Using district fixed effects   ***
****************************************


set scheme plotplainblind

foreach macro in govperform economic dem_support zanu polactivity community state_capacity {

	capture mkdir "plots/district/"
	capture mkdir "plots/district/`macro'"
	
	foreach var of varlist $`macro' {
	
	est clear
	loc j=0
	
	
	
	disp "`var'"
	
		forvalues i = 5(0.5)9.5 {
		
		disp "band = `i'"
		quietly {

			eststo: areg `var' c.treatment native dist_to_border pv ///
			treatment_int native_int pv_int i.round_group year_birth ///
			if dist_to_moz < 150 & abs(dist_to_border) < `i', ///
			absorb(district) cluster(local_auth)
			
			}
			
			est sto m`j'
			loc `++j'
			
			*tab native pv if dist_to_moz < 150 & abs(dist_to_border) < `i'
			*noisily lincom treatment
		}

		coefplot (m0 m1 m2 m3 m4 m5 m7 m6 m7 m8 m9), pstyle(gs15) mcol(black) ///
		keep(treatment) asequation ///
		yline(0, lcolor(black) lwidth(0.25)) levels(90 95) msize(2) symbol(circle) ///
		ciopts(lwidth(0.8 0.5)) yscale(range(-5, 5)) ///
		ylabel(-5(1)5,labcolor(black) labsize(large)) ///
		xlabel(,labcolor(black) labsize(large)) ///
		ytitle("Effect size", axis(1) color(black) size(large)) ///
		xtitle("Bandwidth (km)",  color(black) size(large)) legend(off) ///
		coeflabels(, wrap(15)) ///
		swapnames eqrename (m0 = "5" m1 = "5.5" m2 = "6" m3 = "6.5" m4 = "7" m5 = "7.5" ///
		m6 = "8" m7 = "8.5" m8 = "9" m9 = "9.5") vertical

		graph display Graph, ysize(1) xsize(1)
		graph export "plots/district/`macro'/`var'.pdf", replace

	}	
	
}

**********************************
***   Conley standard errors   ***
**********************************

set scheme plotplainblind

preserve
drop if ward_centroid_y == .
 
foreach macro in govperform economic dem_support zanu polactivity community state_capacity {

	capture mkdir "plots/conley/"
	capture mkdir "plots/conley/`macro'"
	
	foreach var of varlist $`macro' {
	
	est clear
	loc j=0
	
	disp "`var'"
	
		forvalues i = 5(0.5)9.5 {
		
		disp "band = `i'"
		quietly {
	
			eststo: acreg `var' c.treatment native dist_to_border pv ///
			treatment_int native_int pv_int year_birth ///
			if dist_to_moz < 150 & abs(dist_to_border) < `i', ///
			spatial latitude(ward_centroid_x) longitude(ward_centroid_y) ///
			dist(`i') pfe1(provincepc) pfe2(round_group)
			
			}
			
			est sto m`j'
			loc `++j'
			
			*tab native pv if dist_to_moz < 150 & abs(dist_to_border) < `i'
			*noisily lincom treatment
		}

		coefplot (m0 m1 m2 m3 m4 m5 m7 m6 m7 m8 m9), pstyle(gs15) mcol(black) ///
		keep(treatment) asequation ///
		yline(0, lcolor(black) lwidth(0.25)) levels(90 95) msize(2) symbol(circle) ///
		ciopts(lwidth(0.8 0.5)) yscale(range(-5, 5)) ///
		ylabel(-5(1)5,labcolor(black) labsize(large)) ///
		xlabel(,labcolor(black) labsize(large)) ///
		ytitle("Effect size", axis(1) color(black) size(large)) ///
		xtitle("Bandwidth (km)",  color(black) size(large)) legend(off) ///
		coeflabels(, wrap(15)) ///
		swapnames eqrename (m0 = "5" m1 = "5.5" m2 = "6" m3 = "6.5" m4 = "7" m5 = "7.5" ///
		m6 = "8" m7 = "8.5" m8 = "9" m9 = "9.5") vertical

		quietly graph display Graph, ysize(1) xsize(1)
		graph export "plots/conley/`macro'/`var'.pdf", replace

	}	
	
}

restore


*****************************
***   Still rural today   ***
*****************************


set scheme plotplainblind

foreach macro in govperform economic dem_support zanu polactivity community state_capacity {

	capture mkdir "plots/rural/"
	capture mkdir "plots/rural/`macro'"
	
	foreach var of varlist $`macro' {
	
	est clear
	loc j=0
	
	disp "`var'"
	
		forvalues i = 5(0.5)9.5 {
		
		disp "band = `i'"
		quietly {

			eststo: areg `var' c.treatment native dist_to_border pv ///
			treatment_int native_int pv_int i.round_group year_birth ///
			if dist_to_moz < 150 & abs(dist_to_border) < `i' & urban_rural == 2, ///
			absorb(province) cluster(local_auth)
			
			}
			
			est sto m`j'
			loc `++j'
			
			*tab native pv if dist_to_moz < 150 & abs(dist_to_border) < `i'
			*noisily lincom treatment
		}

		coefplot (m0 m1 m2 m3 m4 m5 m7 m6 m7 m8 m9), pstyle(gs15) mcol(black) ///
		keep(treatment) asequation ///
		yline(0, lcolor(black) lwidth(0.25)) levels(90 95) msize(2) symbol(circle) ///
		ciopts(lwidth(0.8 0.5)) yscale(range(-5, 5)) ///
		ylabel(-5(1)5,labcolor(black) labsize(large)) ///
		xlabel(,labcolor(black) labsize(large)) ///
		ytitle("Effect size", axis(1) color(black) size(large)) ///
		xtitle("Bandwidth (km)",  color(black) size(large)) legend(off) ///
		coeflabels(, wrap(15)) ///
		swapnames eqrename (m0 = "5" m1 = "5.5" m2 = "6" m3 = "6.5" m4 = "7" m5 = "7.5" ///
		m6 = "8" m7 = "8.5" m8 = "9" m9 = "9.5") vertical

		quietly graph display Graph, ysize(1) xsize(1)
		graph export "plots/rural/`macro'/`var'.pdf", replace

	}	
	
}


**************************
***   Cohort effects   ***
**************************


set scheme plotplainblind
capture mkdir "plots/age_dist"

forvalues i = 5(0.5)9.5 {
		
	disp "band = `i'"
	
	hist year_birth if dist_to_moz < 150 & abs(dist_to_border) < `i'
	graph display Graph, ysize(1) xsize(1)
	graph export "plots/age_dist/age `i' band.pdf", replace

}



foreach macro in govperform economic dem_support zanu polactivity community state_capacity {

	capture mkdir "plots/`macro'"
	
	foreach var of varlist $`macro' {
	
	est clear
	loc j=0
	
	disp "`var'"
	
		forvalues i = 5(0.5)9.5 {
		
		disp "band = `i'"
		quietly {
			
			ren treatment treat_`j'
			eststo: areg `var' c.treat_`j' native dist_to_border pv ///
			treatment_int native_int pv_int ///
			i.round_group year_birth ///
			if dist_to_moz < 150 & abs(dist_to_border) < `i' & post == 1, ///
			absorb(province) cluster(local_auth)

			est sto m`j'_postwar
			
			eststo: areg `var' c.treat_`j' native dist_to_border pv ///
			treatment_int native_int pv_int ///
			i.round_group male year_birth ///
			if dist_to_moz < 150 & abs(dist_to_border) < `i' & post == 0, ///
			absorb(province) cluster(local_auth)
			
			est sto m`j'_war
			
			ren treat_`j' treatment
			
			}
			
			loc `++j'
			
		}
		
		coefplot (m0_war, offset(-0.075) pstyle(p7) label("Pre-independence") mcol(sea) keep(treat_0)) ///
		(m0_postwar, offset(0.075) pstyle(p3) label("Post-independence") mcol(sky) keep(treat_0)) ///
		(m1_war, offset(-0.075) pstyle(p7) label("Pre-independence") mcol(sea) keep(treat_1)) ///
		(m1_postwar, offset(0.075) pstyle(p3) label("Post-independence") mcol(sky) keep(treat_1)) ///
		(m2_war, offset(-0.075) pstyle(p7) label("Pre-independence") mcol(sea) keep(treat_2)) ///
		(m2_postwar, offset(0.075) pstyle(p3) label("Post-independence") mcol(sky) keep(treat_2)) ///
		(m3_war, offset(-0.075) pstyle(p7) label("Pre-independence") mcol(sea) keep(treat_3)) ///
		(m3_postwar, offset(0.075) pstyle(p3) label("Post-independence") mcol(sky) keep(treat_3)) ///
		(m4_war, offset(-0.075) pstyle(p7) label("Pre-independence") mcol(sea) keep(treat_4)) ///
		(m4_postwar, offset(0.075) pstyle(p3) label("Post-independence") mcol(sky) keep(treat_4)) ///
		(m5_war, offset(-0.075) pstyle(p7) label("Pre-independence") mcol(sea) keep(treat_5)) ///
		(m5_postwar, offset(0.075) pstyle(p3) label("Post-independence") mcol(sky) keep(treat_5)) ///
		(m6_war, offset(-0.075) pstyle(p7) label("Pre-independence") mcol(sea) keep(treat_6)) ///
		(m6_postwar, offset(0.075) pstyle(p3) label("Post-independence") mcol(sky) keep(treat_6)) ///
		(m7_war, offset(-0.075) pstyle(p7) label("Pre-independence") mcol(sea) keep(treat_7)) ///
		(m7_postwar, offset(0.075) pstyle(p3) label("Post-independence") mcol(sky) keep(treat_7)) ///
		(m8_war, offset(-0.075) pstyle(p7) label("Pre-independence") mcol(sea) keep(treat_8)) ///
		(m8_postwar, offset(0.075) pstyle(p3) label("Post-independence") mcol(sky) keep(treat_8)) ///
		(m9_war, offset(-0.075) pstyle(p7) label("Pre-independence") mcol(sea) keep(treat_9)) ///
		(m9_postwar, offset(0.075) pstyle(p3) label("Post-independence") mcol(sky) keep(treat_9)), ///
		yline(0, lcolor(black) lwidth(0.25)) levels(90 95) msize(1.5) symbol(circle) ///
		ciopts(lwidth(0.35 0.2)) yscale(range(-5, 5)) ///
		ylabel(-5(1)5,labcolor(black) labsize(large)) ///
		xlabel(,labcolor(black) labsize(large)) ///
		ytitle("Effect size", axis(1) color(black) size(medlarge)) ///
		xtitle("Bandwidth (km)",  color(black) size(medlarge)) legend(off) ///
		coeflabels(, wrap(15)) ///
		rename (treat_0 = "5" treat_1 = "5.5" treat_2 = "6" treat_3 = "6.5" treat_4 = "7" ///
			treat_5 = "7.5" treat_6 = "8" treat_7 = "8.5" treat_8 = "9" treat_9 = "9.5") vertical

		graph display Graph, ysize(1) xsize(1)
		graph export "plots/`macro'/`var'_post.pdf", replace

	}	
	
}


***********************************
***   With geography controls   ***
***********************************


set scheme plotplainblind

foreach macro in govperform economic dem_support zanu polactivity community state_capacity {

	capture mkdir "plots/geo_controls/"
	capture mkdir "plots/geo_controls/`macro'"
	
	foreach var of varlist $`macro' {
	
	est clear
	loc j=0
	
	
	
	disp "`var'"
	
		forvalues i = 5(0.5)9.5 {
		
		disp "band = `i'"
		quietly {

			eststo: areg `var' c.treatment native dist_to_border pv ///
			treatment_int native_int pv_int i.round_group year_birth ///
			share_crop share_forest share_natpark ///
			if dist_to_moz < 150 & abs(dist_to_border) < `i', ///
			absorb(province) cluster(local_auth)
			
			}
			
			est sto m`j'
			loc `++j'
			
			*tab native pv if dist_to_moz < 150 & abs(dist_to_border) < `i'
			*noisily lincom treatment
		}

		coefplot (m0 m1 m2 m3 m4 m5 m7 m6 m7 m8 m9), pstyle(gs15) mcol(black) ///
		keep(treatment) asequation ///
		yline(0, lcolor(black) lwidth(0.25)) levels(90 95) msize(2) symbol(circle) ///
		ciopts(lwidth(0.8 0.5)) yscale(range(-5, 5)) ///
		ylabel(-5(1)5,labcolor(black) labsize(large)) ///
		xlabel(,labcolor(black) labsize(large)) ///
		ytitle("Effect size", axis(1) color(black) size(large)) ///
		xtitle("Bandwidth (km)",  color(black) size(large)) legend(off) ///
		coeflabels(, wrap(15)) ///
		swapnames eqrename (m0 = "5" m1 = "5.5" m2 = "6" m3 = "6.5" m4 = "7" m5 = "7.5" ///
		m6 = "8" m7 = "8.5" m8 = "9" m9 = "9.5") vertical

		graph display Graph, ysize(1) xsize(1)
		graph export "plots/geo_controls/`macro'/`var'.pdf", replace

	}	
	
}


*******************************************
***   Ward population characteristics   ***
*******************************************

clear
use "data/zim_geo_wards.dta"

gen pv = 0
replace pv = 1 if pv_area != ""
replace pv = 1 if pv_border == 1
encode pv_area, gen(pv_code)

gen rural = strpos(local_auth, "RDC") > 0
replace rural = . if local_auth == ""

gen dist_to_border = native_dis
replace dist_to_border = -native_dis if native == 0
gen treatment = (pv == 1 & dist_to_border > 0) if pv != . & dist_to_border != .
gen treatment_int = treatment * dist_to_border
gen native_int = native * dist_to_border
gen pv_int = pv * dist_to_border
label var treatment "Effect of PV"

global population pct_migran pop popdens12

foreach macro in population {

	capture mkdir "plots/`macro'"
	
	foreach var of varlist $`macro' {
	
	est clear
	loc j=0
	
	disp "`var'"
	
		forvalues i = 5(0.5)9.5 {
		
		disp "band = `i'"
		quietly {

			eststo: areg `var' c.treatment native dist_to_border pv ///
			treatment_int native_int pv_int ///
			if dist_to_moz < 150 & abs(dist_to_border) < `i', ///
			absorb(province) cluster(local_auth)
			
			}
			
			est sto m`j'
			loc `++j'
			
		}

		coefplot (m0 m1 m2 m3 m4 m5 m7 m6 m7 m8 m9), pstyle(gs15) mcol(black) ///
		keep(treatment) asequation ///
		yline(0, lcolor(black) lwidth(0.25)) levels(90 95) msize(2) symbol(circle) ///
		ciopts(lwidth(0.8 0.5)) yscale(range(-5, 5)) ///
		ylabel(-5(1)5,labcolor(black) labsize(large)) ///
		xlabel(,labcolor(black) labsize(large)) ///
		ytitle("Effect size", axis(1) color(black) size(large)) ///
		xtitle("Bandwidth (km)",  color(black) size(large)) legend(off) ///
		coeflabels(, wrap(15)) ///
		swapnames eqrename (m0 = "5" m1 = "5.5" m2 = "6" m3 = "6.5" m4 = "7" m5 = "7.5" ///
		m6 = "8" m7 = "8.5" m8 = "9" m9 = "9.5") vertical

		graph display Graph, ysize(1) xsize(1)
		graph export "plots/`macro'/wardlevel_`var'.pdf", replace

	}	
	
}

